library(png)
format(Sys.time(), format="%d_%b_%Y")
[1] "27_Feb_2017"
## First consolidate the available files into a single table
path <- "~/Box Sync/Four model compare/Module 2 extinct"
setwd(path)
Error in setwd(path) : cannot change working directory
setwd("~/Box Sync/colliding ranges/Simulations_humans/Results/available daily summaries")
The working directory was changed to /Users/Ty/Box Sync/colliding ranges/Simulations_humans/Results/available daily summaries inside a notebook chunk. The working directory will be reset when the chunk is finished running. Use the knitr root.dir option in the setup chunk to change the the working directory for notebook chunks.
details <- file.info(list.files())
trimmed_details <- details[which(list.files() == list.files(pattern = "Four_model_compare_results_extinct*")),]
ord <- order(trimmed_details$mtime, decreasing = TRUE)
rownames(trimmed_details[ord,])[1]
[1] "Four_model_compare_results_extinct_27_Feb_2017.Rdata"
load(rownames(trimmed_details[ord,])[1])
extinct <- Concatenated_data
trimmed_details <- details[which(list.files() != list.files(pattern = "Four_model_compare_results_extinct*")),]
ord <- order(trimmed_details$mtime, decreasing = TRUE)
rownames(trimmed_details[ord,])[1]
[1] "Four_model_compare_results_27_Feb_2017_crop_to_2861.Rdata"
load(rownames(trimmed_details[ord,])[1])
extant <- Concatenated_data
for(i in c(9,10,11,12,14,15,16,17,19,20,21,22,24,25,26,27,29,30,31,32)){
extinct[which(is.nan(as.numeric(as.character(extinct[, i]))) == TRUE), i] <- NA
}
for(i in c(9,10,11,12,14,15,16,17,19,20,21,22,24,25,26,27,29,30,31,32)){
extant[which(is.nan(as.numeric(as.character(extant[, i]))) == TRUE), i] <- NA
}
i <- 19
for(i in c(20,21,24,25,26,27)){
extinct[which(as.numeric(as.character(extinct[, i])) == 0), i] <- NA
}
for(i in c(20,21,24,25,26,27)){
extant[which(as.numeric(as.character(extant[, i])) == 0), i] <- NA
}
xlimit <- c(0,1)
ylimit <- c(0,600)
maincex <- 0.9
png(file="Global_success_rate_per_parameter.png", width=8.5, height=11, units="in", res=300)
par(mfrow=c(5,4), mar=c(3,3,3,0))
hist(as.numeric(as.character(extinct[,9])), main="speciation of F in F env", col=adjustcolor("firebrick", alpha=0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex)
hist(as.numeric(as.character(extant[,9])), main="speciation of F in F env", col=adjustcolor("cornflowerblue", alpha=0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex, add=TRUE)
hist(as.numeric(as.character(extinct[,10])), main="speciation of D in F env", col=adjustcolor("firebrick", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex)
hist(as.numeric(as.character(extant[,10])), main="speciation of D in F env", col=adjustcolor("cornflowerblue", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex, add=TRUE)
hist(as.numeric(as.character(extinct[,11])), main="speciation of F in D env", col=adjustcolor("firebrick", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex)
hist(as.numeric(as.character(extant[,11])), main="speciation of F in D env", col=adjustcolor("cornflowerblue", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex, add=TRUE)
hist(as.numeric(as.character(extinct[,12])), main="speciation of D in D env", col=adjustcolor("firebrick", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex)
hist(as.numeric(as.character(extant[,12])), main="speciation of D in D env", col=adjustcolor("cornflowerblue", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex, add=TRUE)
#######
hist(as.numeric(as.character(extinct[, 14])), main="extinction of F in F env", col=adjustcolor("firebrick", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex)
hist(as.numeric(as.character(extant[, 14])), main="extinction of F in F env", col=adjustcolor("cornflowerblue", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex, add=TRUE)
hist(as.numeric(as.character(extinct[, 15])), main="extinction of D in F env", col=adjustcolor("firebrick", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex)
hist(as.numeric(as.character(extant[, 15])), main="extinction of D in F env", col=adjustcolor("cornflowerblue", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex, add=TRUE)
hist(as.numeric(as.character(extinct[, 16])), main="extinction of F in D env", col=adjustcolor("firebrick", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex)
hist(as.numeric(as.character(extant[, 16])), main="extinction of F in D env", col=adjustcolor("cornflowerblue", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex, add=TRUE)
hist(as.numeric(as.character(extinct[, 17])), main="extinction of D in D env", col=adjustcolor("firebrick", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex)
hist(as.numeric(as.character(extant[, 17])), main="extinction of D in D env", col=adjustcolor("cornflowerblue", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex, add=TRUE)
######
hist(as.numeric(as.character(extinct[, 29])), main="arisal of F in F env", col=adjustcolor("firebrick", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex)
hist(as.numeric(as.character(extant[, 29])), main="arisal of F in F env", col=adjustcolor("cornflowerblue", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex, add=TRUE)
hist(as.numeric(as.character(extinct[, 30])), main="arisal of D in F env", col=adjustcolor("firebrick", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex)
hist(as.numeric(as.character(extant[, 30])), main="arisal of D in F env", col=adjustcolor("cornflowerblue", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex, add=TRUE)
hist(as.numeric(as.character(extinct[, 31])), main="arisal of F in D env", col=adjustcolor("firebrick", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex)
hist(as.numeric(as.character(extant[, 31])), main="arisal of F in D env", col=adjustcolor("cornflowerblue", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex, add=TRUE)
hist(as.numeric(as.character(extinct[, 32])), main="arisal of D in D env", col=adjustcolor("firebrick", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex)
hist(as.numeric(as.character(extant[, 32])), main="arisal of D in D env", col=adjustcolor("cornflowerblue", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex, add=TRUE)
######
hist(as.numeric(as.character(extinct[, 19])), main="NOPE -- Diffusion: source F, target F", col=adjustcolor("firebrick", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= c(0,18000), cex.main= maincex)
hist(as.numeric(as.character(extant[, 19])), main="NOPE -- Diffusion: source F, target F", col=adjustcolor("cornflowerblue", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= c(0,18000), cex.main= maincex, add=TRUE)
hist(as.numeric(as.character(extinct[, 20])), main="Diffusion: source D, target F", col=adjustcolor("firebrick", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex)
hist(as.numeric(as.character(extant[, 20])), main="Diffusion: source D, target F", col=adjustcolor("cornflowerblue", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex, add=TRUE)
hist(as.numeric(as.character(extinct[, 21])), main="Diffusion: source F, target D", col=adjustcolor("firebrick", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex)
hist(as.numeric(as.character(extant[, 21])), main="Diffusion: source F, target D", col=adjustcolor("cornflowerblue", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex, add=TRUE)
hist(as.numeric(as.character(extinct[, 22])), main="NOPE -- Diffusion: source D, target D", col=adjustcolor("firebrick", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= c(0,18000), cex.main= maincex)
hist(as.numeric(as.character(extant[, 22])), main="NOPE -- Diffusion: source D, target D", col=adjustcolor("cornflowerblue", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= c(0,18000), cex.main= maincex, add=TRUE)
####
hist(as.numeric(as.character(extinct[, 24])), main="Takeover: source F, target F", col=adjustcolor("firebrick", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex)
hist(as.numeric(as.character(extant[, 24])), main="Takeover: source F, target F", col=adjustcolor("cornflowerblue", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex, add=TRUE)
hist(as.numeric(as.character(extinct[, 25])), main="Takeover: source D, target F", col=adjustcolor("firebrick", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex)
hist(as.numeric(as.character(extant[, 25])), main="Takeover: source D, target F", col=adjustcolor("cornflowerblue", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex, add=TRUE)
hist(as.numeric(as.character(extinct[, 26])), main="Takeover: source F, target D", col=adjustcolor("firebrick", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex)
hist(as.numeric(as.character(extant[, 26])), main="Takeover: source F, target D", col=adjustcolor("cornflowerblue", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex, add=TRUE)
hist(as.numeric(as.character(extinct[, 27])), main="Takeover: source D, target D", col=adjustcolor("firebrick", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex)
hist(as.numeric(as.character(extant[, 27])), main="Takeover: source D, target D", col=adjustcolor("cornflowerblue", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex, add=TRUE)
dev.off()
null device
1
png(file="extiction minus extant per outcome.png", width=8.5, height=11, units="in", res=300)
par(mfrow=c(3,1))
plot(as.numeric(as.character(extinct[,9])), as.numeric(as.character(extinct[,14])), xlab="speciation", ylab="extinction", col= adjustcolor("firebrick", alpha=0.2), pch=19, cex=0.6, ylim=c(0,1))
plot(as.numeric(as.character(extant[,9])), as.numeric(as.character(extant[,14])), xlab="speciation", ylab="extinction", col= adjustcolor("cornflowerblue", alpha=0.2), pch=19, cex=0.6, ylim=c(0,1))
plot(as.numeric(as.character(extinct[,9])), as.numeric(as.character(extinct[,14])), xlab="speciation", ylab="extinction", col= adjustcolor("firebrick", alpha=0.2), pch=19, cex=0.6, ylim=c(0,1))
points(as.numeric(as.character(extant[,9])), as.numeric(as.character(extant[,14])), xlab="speciation", ylab="extinction", col= adjustcolor("cornflowerblue", alpha=0.2), pch=19, cex=0.6)
dev.off()
null device
1
load('~/Box Sync/colliding ranges/Simulations_humans/Available trees/real.analysis.RData')
#Concatenated_data <- Concatenated_data[Concatenated_data[, 2] == "stats.no.bTO", ]
#Concatenated_data <- Concatenated_data[Concatenated_data[, 6] != "05", ]
# Concatenated_data[, 6] <- as.numeric(Concatenated_data[, 6])
# # Concatenated_data[original[, 2] == "background_takeover", 6] <- Concatenated_data[original[, 2] == "background_takeover", 6] + 4
Concatenated_data[, 6] <- factor(Concatenated_data[, 6])
#head(Concatenated_data)
#names(Concatenated_data)
PCAdata <- Concatenated_data[, -(1:35)]
PCAdata <- PCAdata[, -12]
PCAdata <- apply(PCAdata, 2, as.numeric)
remove <- apply(is.na(PCAdata), 1, any)
PCAdata <- PCAdata[!remove, ]
# Predictions
library(randomForest)
data.analysis.comp2 <- data.frame("Model" = as.factor(Concatenated_data[!remove, 6]),
PCAdata)
data.analysis.comp2$sprate <- data.analysis.comp2$trait_1_speciation/data.analysis.comp2$trait_2_speciation
data.analysis.comp2$extrate <- data.analysis.comp2$trait_1_extinction/data.analysis.comp2$trait_2_extinction
#load("Real_phy/real.analysis.RData")
a <- as.data.frame(real.analysis$results_summary_of_single_value_outputs)
a$sprate <- a$trait_1_speciation / a$trait_2_speciation
a$extrate <- a$trait_1_extinction / a$trait_2_extinction
data.analysis.comp3 <- data.analysis.comp2[, -c(2, 13:14, 16:20)]
#data.analysis.comp3 <- data.analysis.comp3[data.analysis.comp3$Model %in% 1:4, ]
#data.analysis.comp3$Model <- factor(data.analysis.comp3$Model)
#sub <- unlist(lapply(as.list(c(1:4)), function(x, y) {
# sample(which(y$Model == x), min(table(data.analysis.comp3$Model)))},
# y = data.analysis.comp3))
# data.analysis.comp3 <- data.analysis.comp3[sub, ]
fun <- function(x, y, per = .33) {sample(which(y$Model == x), round(table(y$Model)[1]*per))}
sub.test <- unlist(lapply(as.list(paste0(0, c(1:4))), fun,
y = data.analysis.comp3))
test2 <- data.analysis.comp3[sub.test, 2:ncol(data.analysis.comp3)]
test1 <- data.analysis.comp3[sub.test, 1]
train <- data.analysis.comp3[-sub.test, ]
(fit <- randomForest(Model ~ ., data=train, xtest = test2, ytest = test1,
importance=TRUE, ntree=2400, keep.forest = TRUE, replace=TRUE, norm.votes=FALSE))
Call:
randomForest(formula = Model ~ ., data = train, xtest = test2, ytest = test1, importance = TRUE, ntree = 2400, keep.forest = TRUE, replace = TRUE, norm.votes = FALSE)
Type of random forest: classification
Number of trees: 2400
No. of variables tried at each split: 4
OOB estimate of error rate: 36.25%
Confusion matrix:
01 02 03 04 class.error
01 2796 579 473 8 0.2748963
02 669 2642 214 29 0.2566123
03 809 295 1060 25 0.5157606
04 63 441 121 55 0.9191176
Test set error rate: 48.7%
Confusion matrix:
01 02 03 04 class.error
01 1430 258 208 3 0.2469721
02 332 1440 105 22 0.2417062
03 679 258 939 23 0.5055292
04 183 1257 371 88 0.9536598
predictions <- predict(fit,
a,
type="prob")
plot(fit, ylim=c(0,1))

labs <- c("Basic", "+Diffusion", "+Takeover", "+Diffusion +Takeover")
# bar plot
png("Prob_aus.png", width = 25, height = 25, res = 300, units = "in")
par(mar = c(8, 8, 1, 1))
pred <- setNames(as.numeric(predictions), labs)
cols <- rev(c("darkgreen", "red", "blue", "darkorange1"))
barplot(pred, col = cols, ylab = "Proability", cex.lab = 3, cex.names = 2)
dev.off()
null device
1
prop
01 02 03 04
01 62.280130 13.723696 21.829763 4.678007
02 13.289902 64.409881 8.627678 27.764277
03 20.065147 4.757548 53.503185 13.912515
04 4.364821 17.108875 16.039375 53.645200
importance(fit)
01 02 03
Pylo_diversity_is_sum_of_BL 36.27731 26.313769 38.900228
average_phylogenetic_diversity_is_mean_of_BL 42.47717 32.482005 38.519904
variance_Pylo_diversity_is_variance_of_BL 42.15824 16.661635 9.266383
F_quadratic_entropy_is_sum_of_PD 24.95204 19.851036 36.048651
Mean_pairwise_distance 31.03797 25.728248 33.157743
variance_pairwise_distance 31.43706 29.867035 34.057818
Evolutionary_distinctiveness_sum 36.73711 26.605041 37.480353
mean_Phylogenetic_isolation 40.52562 33.527197 41.801064
variance_Phylogenetic_isolation 24.57867 8.325148 37.260442
gamma 26.36037 21.591877 63.126688
extinction_per_speciation 7.99549 41.950330 54.475002
transition_from_trait_1_to_2 22.53100 22.856490 41.900317
transition_from_trait_2_to_1 33.37167 33.715083 45.090491
transition_rate_ratio_1to2_over_2to1 50.20133 36.600927 35.430284
Phylogenetic_signal 64.80645 43.118360 53.232738
spatial.tests.fora 43.20987 86.600931 89.746230
spatial.tests.dom 82.03166 30.318014 47.558336
prevalence 54.14112 9.616747 35.927323
sprate 41.34870 28.286586 85.062626
extrate 28.88581 24.078228 21.624649
04 MeanDecreaseAccuracy
Pylo_diversity_is_sum_of_BL 45.76293 62.84990
average_phylogenetic_diversity_is_mean_of_BL 32.73694 67.29494
variance_Pylo_diversity_is_variance_of_BL 27.05678 54.35726
F_quadratic_entropy_is_sum_of_PD 43.18824 68.31908
Mean_pairwise_distance 42.30806 63.22960
variance_pairwise_distance 22.73712 60.82865
Evolutionary_distinctiveness_sum 44.98800 63.01409
mean_Phylogenetic_isolation 33.38172 68.79250
variance_Phylogenetic_isolation 22.99055 49.94780
gamma 51.68400 89.28718
extinction_per_speciation 54.09510 75.50858
transition_from_trait_1_to_2 41.11630 58.91197
transition_from_trait_2_to_1 40.76374 67.48847
transition_rate_ratio_1to2_over_2to1 42.40223 85.56456
Phylogenetic_signal 40.17895 85.09066
spatial.tests.fora 76.44656 125.15481
spatial.tests.dom 53.97271 104.59956
prevalence 58.91447 79.90220
sprate 36.41170 92.93119
extrate 36.77120 58.20976
MeanDecreaseGini
Pylo_diversity_is_sum_of_BL 210.2330
average_phylogenetic_diversity_is_mean_of_BL 182.1725
variance_Pylo_diversity_is_variance_of_BL 164.4793
F_quadratic_entropy_is_sum_of_PD 185.9903
Mean_pairwise_distance 206.3388
variance_pairwise_distance 183.6735
Evolutionary_distinctiveness_sum 207.9225
mean_Phylogenetic_isolation 184.8331
variance_Phylogenetic_isolation 167.5043
gamma 213.3865
extinction_per_speciation 251.9260
transition_from_trait_1_to_2 255.1930
transition_from_trait_2_to_1 262.6217
transition_rate_ratio_1to2_over_2to1 211.2503
Phylogenetic_signal 303.1266
spatial.tests.fora 395.4726
spatial.tests.dom 250.7786
prevalence 226.0943
sprate 244.7882
extrate 190.3293
# Variables importance
imp <- importance(fit)
imp <- apply(imp, 2, function(x) (x - min(x))/(max(x) - min(x)))
imp <- imp[sort(imp[, 5], index.return = TRUE, decreasing = TRUE)$ix, ]
names <- rownames(imp)
names[names == "spatial.tests.fora"] <- "Space F"
names[names == "spatial.tests.dom"] <- "Space D"
names[names == "sprate"] <- "Sp(ratio)"
names[names == "transition_from_trait_1_to_2"] <- "TR(1-2)"
names[names == "transition_from_trait_2_to_1"] <- "TR(2-1)"
names[names == "Phylogenetic_signal"] <- "PhySig(D)"
names[names == "Evolutionary_distinctiveness_sum"] <- "EDsum"
names[names == "Pylo_diversity_is_sum_of_BL"] <- "PDsum"
names[names == "transition_rate_ratio_1to2_over_2to1"] <- "TR(ratio)"
names[names == "gamma"] <- "Gamma"
names[names == "mean_Phylogenetic_isolation"] <- "MPI"
names[names == "extrate"] <- "Ext(ratio)"
names[names == "average_phylogenetic_diversity_is_mean_of_BL"] <- "PDmean"
names[names == "extinction_per_speciation"] <- "DR"
names[names == "variance_Phylogenetic_isolation"] <- "VPI"
names[names == "F_quadratic_entropy_is_sum_of_PD"] <- "F"
names[names == "Mean_pairwise_distance"] <- "MPD"
names[names == "variance_Pylo_diversity_is_variance_of_BL"] <- "PDvar"
names[names == "variance_pairwise_distance"] <- "VPD"
png("var_import_all.png", width = 25, height = 25, unit="in", res=300)
par(mar = c(10, 18, 1, 1))
plot(x = rev(imp[, 5]), y = 1:nrow(imp), type = "l", yaxt = "n",
ylab = "", xlab = "Variable Importance",
xlim = c(0, 1), lwd = 2, cex.lab = 4)
for (i in 1:nrow(imp)) {
abline(h = i, lty = 3, col = "gray80")
}
abline(v = seq(0, 1, 1/19), lty = 3, col = "gray80")
lines(x = rev(imp[, 4]), y = 1:nrow(imp), col = "darkgreen", lwd = 2)
lines(x = rev(imp[, 3]), y = 1:nrow(imp), col = "red", lwd = 2)
lines(x = rev(imp[, 2]), y = 1:nrow(imp), col = "blue", lwd = 2)
lines(x = rev(imp[, 1]), y = 1:nrow(imp), col = "darkorange1", lwd = 2)
lines(x = rev(imp[, 5]), y = 1:nrow(imp), lwd = 3)
points(x = rev(imp[, 4]), y = 1:nrow(imp), col = "darkgreen", cex = 2)
points(x = rev(imp[, 3]), y = 1:nrow(imp), col = "red", cex = 2)
points(x = rev(imp[, 2]), y = 1:nrow(imp), col = "blue", cex = 2)
points(x = rev(imp[, 1]), y = 1:nrow(imp), col = "darkorange1", cex = 2)
points(x = rev(imp[, 5]), y = 1:nrow(imp), pch = 20, cex = 3)
text(y = 1:nrow(imp), x = par("usr")[1] - .17, labels = rev(names),
srt = 0, pos = 4, xpd = T, cex = 4)
dev.off()
null device
1